rm(list=ls(all=TRUE))

library(MASS)
library(mc2d)

## prob. of 1 for randomizing device or group status indicator ##

p <- 0.25

n.MC <- 2000

pars.MC <- matrix(NA,nrow=n.MC,ncol=3)
mu.hat <- rep(NA,length.out=n.MC) 
delta  <- rep(NA,length.out=n.MC)

n.sim <- 5000  ## sample size for each Monte Carlo simulation ##

X.mu  <- c(.5, 1.5)
X.sig <- rbind(c(.1,.1),c(.1,.2))

Beta.act  <- c(-1,2,-1.3)

for (k in 1:n.MC) {

X<-mvrnorm(n.sim,X.mu,X.sig)
X<-cbind(1,X)

mu.act <- (1+exp(-X%*%Beta.act))^-1
theta  <- rbinom(n=n.sim,size=1,mu.act)

Mprob <- p*theta+(1-p)*(1-theta)

Y<-rbinom(n=n.sim,size=1,Mprob)

## Set up function to calculate parameter values according to E-M algorithm ##

EM <- function(par) {
	
Beta  <-par[1:length(par)]
mu    <- (1+exp(-X%*%Beta))^-1 

Pr.2.3 <- as.numeric(Y==0)*( ((1-p)*mu)/ ((1-p)*mu + p*(1-mu)) ) + as.numeric(Y==1)*( (p*mu)/(p*mu+(1-p)*(1-mu)) )

V <- mu*(1-mu)
X.til<-matrix(NA,nrow=nrow(X),ncol=ncol(X))
for (i in 1:n.sim) {X.til[i,]<-V[i]*X[i,]} 

D <- Pr.2.3 - mu

Beta.j <- Beta + solve(t(X)%*%X.til,tol=.Machine$double.eps^2)%*%t(X)%*%D
Beta.j}

tol <- .Machine$double.eps^.5 # tolerance for stopping algorithm

M<-4000
K<-ncol(X)
init<-rep(0,K) # initial parameter estimates
pars.old <- EM(init)


for (j in 1:M){

pars.new <- EM(pars.old)

abs.diff <- abs(pars.new-pars.old)

if (max(abs.diff) < tol ) break

pars.old <- pars.new}

pars.MC[k,] <- pars.new

X2.l <- quantile(X[,3],.05)
X2.h <- quantile(X[,3],.95)

X.1 <- cbind(X[,1:2],X2.h)
X.0 <- cbind(X[,1:2],X2.l)

mu.ca   <- (1+exp(-X%*%pars.new))^-1
mu.1    <- (1+exp(-X.1%*%pars.new))^-1
mu.0    <- (1+exp(-X.0%*%pars.new))^-1 

mu.hat[k]<-mean(mu.ca)
delta[k]<-mean(mu.1-mu.0)

}

beta0.avg<-mean(pars.MC[,1])
beta1.avg<-mean(pars.MC[,2])
beta2.avg<-mean(pars.MC[,3])

beta0.bias <-beta0.avg-Beta.act[1]
beta1.bias <-beta1.avg-Beta.act[2]
beta2.bias <-beta2.avg-Beta.act[3]

beta0.mse <-mean((pars.MC[,1]-Beta.act[1])^2)
beta1.mse <-mean((pars.MC[,2]-Beta.act[2])^2)
beta2.mse <-mean((pars.MC[,3]-Beta.act[3])^2)

BETAS<-c(beta0.avg,beta1.avg,beta2.avg)
BIAS<-c(beta0.bias,beta1.bias,beta2.bias)
MSE<-c(beta0.mse,beta1.mse,beta2.mse)

QOIs<-cbind(mu.hat,delta)


save(QOIs,file="SSTQOIs5000.RData")

